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AbstractHWe  derive  and  analyze  several  methods  for  systems  of  hyperbolic  equations  with 
wide  ranges  of  signal  speeds.  These  techniques  are  also  useful  for  problems  whose  coefficients  have 
large  mean  values  about  which  they  oscillate  with  small  amplitude.  Our  methods  are  based  on 
additive  splittings  of  the  operators  into  components  that  can  be  approximated  independently  on  the 
different  time  scales,  some  of  which  arc  sometimes  treated  exactly.  The  efficiency  of  the  splitting 
methods  is  seen  to  depend  on  the  error  incurred  in  splitting  the  exact  solution  operator.  This  is 
analyzed  and  a  technique  is  discussed  for  reducing  this  error 'through  a  simple  change  of  variables. 
A  procedure  for  generating  the  appropriate  boundary  data  for  the  intermediate  solutions  is  also 
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1.  Introduction. 


Splitting  methods  lor  time-dependent  partial  differential  equations  have  been  niost  Frequently 
studied  in  the  context  of  spatial  splittings,  as  in  the  approximate  factorization  techniques  for 
ellieiontly  implementing  implicit  algorithms  in  more  than  one  space  dimension[fij,  [llj,  [13].  Some 
attention  has  also  been  given  to  splitting  or  fractional  step  methods  for  problems  where  the 
differential  operator  is  split  up  into  pieces  corresponding  to  different  physical  processes  which  are 
most  naturally  handled  by  different  techniques.  This  has  been  done,  for  example,  with  eonvection- 
dilTusion  and  the  Navicr-Stokcs  cquations(l|,  [4],  [5[. 

More  generally,  a  splitting  method  may  be  useful  any  time  one  is  faced  with  a  problem 


ut  =  Au 


(1.1) 


where  A  is  some  differential  operator  of  the  form 

A  —  A\  +  A% 


such  that  the  problems 


«e  =  Aiu 


(1.2a) 


and 

ut  =  A^u  (1.26) 

are  each  easier  to  solve  than  the  original  problem.  By  alternating  between  solving  (1.2a)  and  (1.2b) 
we  hope  to  compute  a  satisfactory  solution  to  (l.l). 

In  this  paper  we  consider  such  methods  applied  to  a  one-dimensional  quasilinear  hyperbolic 
•  system 

ut  —  A(x,  t,  u)ux  (1.3) 


where  A  is  an  re  X  «  matrix  with  real  eigenvalues.  Wc  assume  that  A  is  of  the  form 


A  =  Af  +  Aa.  (1.4) 

In  our  notation,  “/”  and  “a"  stand  for  “fast”  and  “slow”  respectively,  which  reflects  a  common 
situation  in  which  the  solution  contains  waves  traveling  at  quite  different  wave  speeds.  If  A  is 
constant  then  the  solution  operator  for  the  problem  (1-3)  on  a  single  Liincstep  of  size  k  is  cxp(kAdx), 
that  is  to  say 

u(x,  t  +  k)  =  exp(kAdx)u(x,  l  (1.5) 

For  nonconstant  A  the  solution  operator  is  more  complicated.  Oe  inalysis  will  be  concerned 
mostly  with  the  constant  coefficient  case,  so  we  will  use  the  notation  of  (1.5)  throughout.  The 
ideas  generalize  easily,  but  are  most  intuitively  seen  in  terms  of  exponentials. 

The  additive  splitting  (1.4)  comes  into  play  when  the  solution  operator  cxp(kAdx)  is  ap¬ 
proximated  by  the  product  of  the  solution  operators  for  the  subproblcms 


ut  =  A/ux 


(1.6a) 


and 

Wo  replace  ( 1 .5)  by 


ut  =  Aaux. 


u(x,  t  +  k)  fsa  expjfcyl/t)*)  expjfcyl,^*)^!,  t). 

1 


(1.66) 


An  approximation  lo  u(x,t  +  k )  is  Unis  obtained  by  first,  solving  (l.(ib)  witii  u(x,t)  as  initial  data, 
and  using  the  resulting  solution  as  initial  data  for  (l. (in).  If  A f  and  A*  commute,  this  splitting 
is  exacL.  When  they  do  not  commute,  vve  have  introduced  an  error  which  is  ()(k~).  As  noted  by 
Strang[l6j,  this  error  can  be  reduced  to  0(k '*)  by  use  of  the  splitting 

exp(fc(/ly  +  Aa)()x)  rs  cxp{^kAf()x)cxp(kAa()x)o.xp(^kAfdT).  (1.7) 

Analogous  results  hold  for  the  corresponding  splittings  with  variable  coefficients.  Computations 
confirm  that  the  global  error  is  also  improved  (from  O(k)  to  ()(k2))  by  the  use  of  this  splitting. 

The  numerical  approximations  to  the  solution  operators  cxp(A:A3<9x)  and  cxp(^kA/Ox)  will  be 
denoted  by  Qa(k)  and  Q/(k/‘l)  respectively.  The  numerical  method  based  on  the  Strang  splitting 
(1.7)  is  then 

=  Qf(k/2)Qa(k)Qf(k/2)U^  (1.8) 

where  f/£,  is  the  numerical  approximation  to  u(mh,nk)  on  a  grid  with  Ax  =  h  and  At  =  k.  When 
splitting  a  multidimensional  problem  into  one-dimensional  subproblcins,  this  sort  of  a  splitting 
gives  rise  to  the  so-cal'cd  locally  one  dimensional  (LOD)  method,  a  spatially  split  scheme.  In  the 
present  context  we  will  refer  to  (1.8)  as  the  time-split  method. 

In  practice  Un+l  can  be  computed  via  the  sequence 

U'm  =  Q/(k/2)U^ 

=  Q*W'm  (1.9) 

tC+1  =  Q/W  2)tC 

although  iL  should  be  noted  that  when  several  steps  of  (1.8)  are  applied  successively  the  adjacent 
Q[{k/ 2)  operators  can  be.  combined  into  Qj{k),  and  the  half-step  operators  need  only  be  applied 
at  the  beginning  and  immediately  before  printout,  i.e. 

-  Qf(k(m,{k)Ql(k)---Qf(k)Qa(k)Qf(k/2)U°m. 

There  arc  several  situations  in  which  the  use  of  the  time  split  method  may  lead  to  a  more 
efficient  solution  of  the  original  problem.  We  will  mention  three  such  cases  here.  Our  analysis  will 
be  mostly  concerned  with  the  first  and  last  of  these. 

Problem  1:  Suppose  the  solution  to  (1.3)  contains  both  fast  waves  and  slow  waves,  i.e.  the 
eigenvalues  Hi  of  A  satisfy 

ImiI  <  M  <  •••  <  l/*p|  <  I/vmI  <  •••  <  M- 

Assume  also  that  there  are  relatively  few  elements  of  A  which  contribute  to  the  fast  waves.  We  can 
take  advantage  of  this  structure  by  splitting  the  operator  into  slow  and  fast  parts  and  using  small 
time  steps  only  on  the  fast  part.  That  is,  we  can  choose  k  so  that  exp(&/t5<9x)  can  be  adcquatly 
represented  by  a  single  step  of  some  finite  difference  scheme  and  then  approximate  cxp(gfc/l/dx) 
by  several  steps  of  a  difference  scheme  with  a  smaller  timestep.  Similarly,  we  can  handle  more 
than  two  clusters  of  wave  speeds  by  means  of  further  splittings. 

Such  a  splitting  method  requires  less  work  than  using  small  timesteps  on  the  full  unsplit 
problem,  and  will  thus  be  more  efficient  provided  the  accuracy  is  not  too  adversely  affected  by  the 
error  in  the  splitting.  We  will  sec  that  this  depends  very  much  on  the  problem  at  hand.  In  cases 
where  the  splitting  error  is  small,  the  time-split  method  actually  may  be  more  accurate,  since  we 


will  1)0  able  to  uso  nearly  optimal  mosli  ratios  for  oaeli  duster  to  minimize  the  trimoation  error 
and  improve  other  characteristics  of  the  method,  such  as  its  dissipative  behavior. 

Similar  splittings  have  been  considered  by  Kngquist,  Custafsson  and  Vroeburg[-l]  lor  this  type 
ol"  problem.  However,  the  splitting  in  their  problem  involved  little  interaction  between  the  different 
time  scales,  so  that  many  of  the  problems  we  shall  encounter  were  not  present. 

Example  l.l.  Consider  a  block  triangular  system  of  the  form 


A  = 


A 12 
Az2 


(1.10) 


with  1 1 ^1 1 1 1 1  ~  H/I22II  ~  1  Jind  1.  It  is  reasonable  to  take 


Af  = 


Mu  0 

,  A,  = 

0 

A12 

0  0 

1  J  *5 

0 

A'i-i 

(1.11) 


For  this  problem  the  effectiveness  of  the  split  method  depends  greatly  on  the  coupling  A\n  between 
the  different,  time  scales.  This  is  analyzed  in  section  2  where  we  also  present  a  simple  procedure 
lor  changing  variables  to  reduce  the  coupling. 

Problem  2s  Consider  the  same  situation  as  in  Problem  1,  but  where  the  fast  waves  are 
known  to  be  absent  from  the  physical  solution  of  interest.  Recently  Kreiss[9],  [10]  and  Browning, 
Kasahara,  and  Kreiss[2]  have  considered  some  new  approaches  for  this  problem  which  rely  upon 
properly  preparing  the  data  so  that  the  fast  wave  components  are  eliminated.  These  can  be 
considered  as  projection  techniques.  Majda[12j  has  considered  using  filters  to  suppress  the  fast 
waves  in  the  same  context. 

In  this  case  the  true  solution  should  satisfy 

cxp(kAi)x)u(x,  t)  =  exp(kAadx)u(x,  t) 


providing  the  splitting  between  fast  and  slow  scales  is  done  correctly.  For  variable  coefficient 
problems  it  will  not  be  possible  to  have  the  correct  splitting  at  all  times  ami  .he  operator  Af 
cannot  be  dropped  entirely.  However,  we  can  consider  using  the  time-split  method  (1.8)  with  a  less 
accurate  scheme  for  Qf{k/ 2)  than  is  used  for  Qa(k)  ,  perhaps  by  using  the  same  timestep  for  both 
with  a  larger  spatial  step  for  Qf(k/2)  .  In  such  a  manner  it  may  again  be  possible  to  obtain  the 
same  accuracy  more  efficiently.  Turkel  and  Zwas[18]  have  considered  a  method  for  this  problem 
which  is  similar  in  spirit. 

Problem  3s  Suppose  that  the  coefficients  in  the  problem  (1.3)  have  large  mean  values 
about  which  they  oscillate  with  small  amplitude.  In  this  case  it  may  be  possible  to  split  out  a 
constant  coefficient  problem  which  can  be  solved  exactly,  leaving  behind  the  small  perturbations 
for  Aa.  Then  (1.8)  can  be  used  with  some  large  timestep  approximation  for  Qa(k)  while  Q f[k/ 2)  = 
cxp(^kA/dx)  exactly.  This  is  clearly  more  efficient  than  using  small  timesteps  on  the  unsplit 
problem.  Moreover,  since  the  dominant  part  of  the  operator  is  being  handled  exactly,  great 
increases  in  accuracy  are  also  possible. 

Example  1.2.  The  simplest  example  is  the  scalar  problem 


«t  =  (1  +  a(z))u«  (1.12) 

where  |a(x)|  1  and  we  use  the  splitting 


Af  =  1,  A,  =  ot(x). 
3 


.JIM  IJ  ^III 


Take  k  =  p/»  Tor  some  integer  />.  Tim  operator  <'xp(  f,  k/\ /i)x)  in  known  exactly:  oxp(  f^kA ft),)n(x,  t)  — 
u( x  +■  r,ph,t).  IT  l.ax-Wcndrolf  is  used  lor  the  remaining  suhprohlem  ut  =  x  then  the  method 

(1.8)  can  be  written  as  a  single  stop  method 


IJn  +  l  =  lJn  4- 
u  m  —  u  m  +  p  r  2 


MtC+P+.  -tC+p-.)  +  ip2a(^m)(M:rr 

-  («(xm)  +  tt(xm  -  fc))(^m  +  p  -  Vm+p-i)) 


+  h)  +  «(im)!(^m  +  p+l 


-  U 


m+p‘ 


where  xm  =  xm  +  £ ph .  Notice  that  even  though  this  is  a  scalar  problem,  the  operators  <)x  and 
o(x)dx  do  not  commute  and  so  the  Strang  spli ting  must  be  used. 

The  shallow  ‘niter  equations  provide  a  more  interesting  example  of  Problem  3.  These  are 
discussed  in  section  8. 


General  considerations.  The  effectiveness  of  the  time-split  method  depends  on  the  error 
in  the  splitting  (1.7).  If  this  is  exact  then  the  splitting  method  is  clearly  more  efficient  for  the 
types  of  problems  we  are  considering.  On  the  other  hand,  if  the  splitting  error  dominates,  it 
may  be  necessary  to  reduce  the  timestep  considerably,  eliminating  the  possible  benefits  of  the 
splitting.  In  the  next  section  we  derive  an  explicit  expression  for  the  splitting  error  and  indicate 
how  to  determine  whether  a  splitting  method  is  useful  on  a  given  problem.  We  will  show  that  the 
equations  sometimes  need  to  be  transformed  to  reduce  the  linkage  between  fast  and  slow  modes  in 
order  to  achieve  the  desired  accuracy. 

The  Q  operators  in  the  time-split  method  can  consist  of  one  or  more  steps  of  any  explicit  or 
implicit  scheme  using  two  time  levels.  It  is  not  immediately  clear  how  a  scheme  using  more  than 
two  time  levels  (such  as  leapfrog)  could  be  used.  For  suppose  we  want  to  use  leapfrog  as  Q3[k)  to  go 
from  U'  to  U"  in  (1.9).  Then  we  would  need  some  approximation  to  exp(—kA3i)z)U’  (which  is  not 
Un)  to  play  the  role  of  U ’  at  the  previous  Lime  level.  As  a  first  step  towards  incorporating  multi¬ 
level  schemes  into  the  splitting  framework,  section  \  introduces  a  different  type  of  split  scheme 
which  does  use  leapfrog  for  Q3(k).  This  method  is  based  upon  approximations  to  the  variation 
of  parameters  formula,  or  Duhaincl’s  principal,  and  will  be  called  the  Leapfrog  Duhamcl  method. 
Similar  ideas  have  been  used  for  ordinary  differential  equations  by  Ccrtainc[3).  The  accuracy  and 
stability  of  the  Leapfrog  Duhamcl  method  is  considered  in  sections  5  and  6. 

The  initial  boundary  value  problem  is  considered  in  section  7.  In  most  cases  boundary  data 
will  have  to  bo  supplied  for  the  intermediate  solutions  U  and  U  in  (1.9).  We  consider  the 
problem  of  approximating  the  correct  boundary  data  in  terms  of  the  given  boundary  conditions. 

Some  further  examples  of  splittings  and  computational  results  arc  presented  in  section  8. 


2.  Accuracy  of  the  time-split  method. 

In  this  section  we  consider  discretizations  of  the  approximate  splitting 

u{x,t  +  k)  exp(  jAA fdx)  cxp(kAadz)  cxp(  g  kA/dx)  u{x,t)  (2.1) 

for  the  solution  of  ut  =  Aux  =  (A/  +  A„)ux  with  ||A,||  ||A/||.  Up  until  section  7  we  will  deal 

only  with  the  Cauchy  problem,  where  — oo  <  x  <  oo.  Of  course  these  results  also  hold  for  a  strip 
problem  with  periodic  boundary  conditions,  e.g.,  0  <  x  <  t  and  u(0,  t)  =  u(l,f).  Wc  will  assume 
that  A/  and  A,  arc  constant  matrices  but  our  approach  carries  over  for  more  general  problems 
if  the  exponentials  in  (2.1)  are  replaced  by  the  appropriate  solution  operators.  For  example,  the 
splitting  error  for  the  problem  (1.12)  is  given  in  example  8.2  of  section  8. 

If  A/  and  A,  commute  then  the  splitting  (2.1)  is  exact.  Otherwise  we  define  the  splitting  error 


\ 


operator  /?spi u(fc)  by 


E spiu(A)  =  exp(.jjA:/l/t)z)oxp(A:/i,ai)cxp(^fcA/c)x)  -  exp (fc(A/  +  /13)<9Z) 

=  -  ±k*(\A‘jAa  -  i AjA.Af  +  \A,Aj  (2-2) 

-  {A2aAf  +  AaAfAa  -  I AfA2)03x  +  G{k*). 

The  local  truncation  error  operators  for  the  approximate  solution  operators  Qf(k/ 2)  and 
Qa(k)  are  defined  by 

Ef{k/ 2)  =  Qf( k/2)  -  exp^A/l/d*) 

E.{k)  =  Q,{k)  -  exp (kA,dx). 

Note  that  for  a  second  order  scheme  such  as  Lax-Wendroff  these  are  0{k3).  We  can  now  compute 
the  truncation  error  operator  for  the  split  scheme.  The  numerical  solution  operator  is 

Qf(k/2)Qs(k)Q{(k/2)  =  (c\p(^kAfOx)  +  Ef(k/2))(c\p{kAa0x)  +  Es(k)) 

X  (c\p(^kAfdx)  +  Ef{k/ 2)) 

=  exp(  5  A/l fdx)  cxp(kAadx)  cxp(  \ kA  fdx) 

+  Ea(k)  +  2Ef{k/2)  +  0(k 4) 

=  cxp(A(yt/  +  Aa)dx)  +  Eapiit{k) 

+  E.(k)  +  2Ef(k/2)  4-  0(fc4). 

The  truncation  error  operator  for  the  time-split  method  (TSM)  is  thus 

ET°M[k)  =  Qf(k/2)Qa(k)Qf(k/2)  -  exp (k(Af  +  Aa)dx) 

=  tf.Piu(*)  +  E,(k)  +  2Ef(k/2)  +  0(k*).  [  ' . ' 

For  a  given  problem  this  error  can  be  computed  directly  and  used  to  assess  the  efficiency  of  the 
time-split  method  relative  to  an  unsplit  method. 

In  order  to  illustrate  some  of  the  properties  of  this  method  and  the  elTect  the  splitting  error 
^spiit(^)  has  on  its  utility,  we  restrict  our  attention  to  the  case  where  Qa(k)  consists  of  a  single 
step  of  Lax-Wendroff.  For  convenience  we  use  LW[A,k)  to  denote  the  Lax-WcndrolT  operator, 

LW[A,  k)  =  /  +  kAD0  +  %k2A2D+D- 


where  Do,  D+  and  />_  are  the  standard  centered,  forward,  and  backward  difference  operators, 
respectively.  We  thus  have 


Q»{k)  —  LW[Aa,k). 

(2.4a) 

For  Q/(k/ 2)  we  consider  both 

Qf(k/ 2)  =  oxp{^kAfdx) 

(2.46) 

and 

Q,(k/2)  =  {LW(Af,k/m)r'2 

(2.4c) 

for  some  even  integer  m.  The  situation  (2.4b)  occurs  when  the  solution  operator  exp(£fc/l/ds)  is 
known  exactly,  as  in  Problem  3.  In  (2.4c)  Q/(k/2)  consists  of  m/2  steps  of  Lax-Wendroff  with 
timestep  k/m.  This  might  be  appropriate  when  solving  Problem  1,  for  example. 
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The  standard  error  analysis  lor  Lax- Wend  roll"  shows  that  for  (2.4a)  wo  have 


Ea(k)  =  ^w(k)  =  -  &k*A*  -  kh*A.)Ol  +  O(k'). 


(2.5a) 


When  (2.4b)  is  used  there  is  no  error  on  the  fast  scale  and 

Ef[k/2)  =  Eejxp(k/2)  =  0.  (2.56) 

Otherwise,  when  (2.4c)  is  used, 

Qf(k/2)  =  ([AV(Aflk/m))m 

=  (cxp(£'V5*)  -  i(&*f  ~  &h*A/)di)  ' 

=  e*j>(£kAfOx)  -  -  lfrAf^dl  +  0(k 4) 

and  so 

EM 2)  =  EfW(k/2)  =  -  M&A*f  -  kh*Af)dl  4-  0(*4).  (2.5c) 

In  this  case  we  must  choose  an  appropriate  value  of  m,  the  number  of  small  timesteps  used  within 
each  large  timestep.  For  fixed  h,  the  error  Ehw( k/2)  does  not  approach  zero  as  m  — *•  oo.  From 
(2.5c)  it  seems  unreasonable  to  take  m  any  larger  than  a  value  for  which  j|  ~iAf\\  «  \\kh2Af\\. 
This  suggests  taking 

m  as  -|(y4/||-  (2.6) 

The  proper  choice  of  m  may  also  be  infiucnccd  by  stability  requirements.  Determining  the  stability 
of  the  operator  Q  f(k/2)Qa(k)Q  f(k/2)  is  in  general  a  difficult  problem,  which  will  be  considered 
in  some  detail  in  section  5.  It  will  be  shown  there  that  for  some  problems  the  product  operator 
is  stable  provided  Qf(k/ 2)  and  Qa(k)  are  each  stable  independently.  It  is  well  known  that  for 
Lax-WcndrolT  the  stability  condition  on  Q/(k/2)  is  p(Af)k/rnh  <  1,  i.e.,  m  >  p(Af)k/h.  The  m 
given  in  (2.6)  is  consistent  with  this  requirement.  Also  note  that  for  k/h  as  l/||/\s||,  (2.6)  becomes 

m  »  WAfWlWMl-  _ 

When  the  splitting  error  Eap\\^(k)  is  negligible  compared  to  the  other  terms  in  (2.3),  the 
truncation  error  for  the  split  scheme  becomes 


ETSM(k)  =  E.(k)  +  2Ef(k/2)  +  0(k4) 

=  -  M(Al  +  £ A ))  -  kh*(Aa  +  Af))d3x  +  0(k4). 


This  error  is  roughly  the  same  as  we  would  obtain  using  (LW(A,  k/m))m,  i.e.,  Lax-Wcndroff  with 
'small  steps  on  the  unsplit  problem.  The  truncation  error  Tor  the  unsplit  method  can  be  derived  in 
the  same  manner  as  (2.5c)  to  obtain 

ELW{k)  =  -  &£{A,  +  Aaf  -  kh*(A,  +  A,))dl  +  0(k 4) 

»  -  K&tf  -  kh2(A/  +  A*))dl 

Thus  we  do  almost  as  well  by  taking  large  steps  with  Aa  and  small  steps  with  Af  as  we  would  by 
taking  small  steps  on  the  unsplit  problem.  This  can  lead  to  considerable  savings.  If  Qf[k/ 2)  = 
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cxp(}jkAfDx)  the  results  are  even  more  striking;.  Now  the-orror  (2..'I)  is  simply 

tiT*“{k)  =  /?a(fc)  +  Oik*)  =  -  &k3s\]  -  kh2A.)Ol  +  0(k*). 


Comparing  this  to  (2.7)  shows  that  the  split  scheme  is  considerably  more  accurate.  It  also  requires 
less  work,  since  now  nothing  is  computed  using  small  steps. 

The  results  of  the  last  paragraph  are  all  based  on  the  assumption  that  /?spiit(fc)  is  negligible. 
In  practice  /i’gI>iit(fc)  may  easily  dominate  the  discretization  error  lCa[k)  +  2ISf(k/2).  in  this  case  the 
split  scheme  is  less  accurate  than  Lax-Wendroff  with  timestep  k/m.  Nonetheless  the  split  scheme 
may  be  preferable.  It  may  be  possible  to  use  the  split  scheme  with  smaller  k  and  li  to  obtain 
better  accuracy  while  still  requiring  less  work  than  the  unsplit  scheme.  The  proper  quantity  for 
comparison  is  the  work  required  to  obtain  a  given  accuracy.  This  can  be  estimated  and  compared 
for  various  methods  as  we  now  do.  Under  some  mild  assumptions,  we  will  see  that  the  methods 
(2.4a,b)  and  (2.1a,c)  arc  always  more  efficient  than  the  unsplit  scheme  (providing  we  choose  k/h 
properly). 

Work  comparisons.  We  will  compute  expressions  Tor  the  work  required  to  obtain  a  solution 
at  time  t  —  1  with  error  at  most  t.  All  of  the  bounds  below  are  rough  order  of  magnitude  bounds 
which  are  sufficient  for  our  present  purpose.  Suppose  that 

.  ||4|  «  \\Af\\  =  a,  ||Aa||  =  b 

where  b/a  —  e  1.  Also  suppose  that  ||us,x)j  **  1.  This  is  for  convenience  only,  since  it  removes 
one  common  factor  from  all  of  the  bounds  below. 

We  will  first  analyze  the  unsplit  Lax-Wendroff  method  LW(A,k).  Suppose  that  W  is  the 
work  required  to  compute  LW[A,  at  a  single  point  xm.  Then  the  work  required  to  advance 

the  solution  on  a  unit  i-interval  by  one  unit  of  time  is  W/kh  =  \W/k2  if  k  =  X/i.  The  error 
committed  in  one  time  unit  using  the  unsplit  method  is  bounded  by 

||((LW(/l,A:))*/fc-exP(/iai))«|! 

<  l(&3M\3  +  kh2\\A\\)  +  0(k<)) 

<  ±k2(a3  +  a/\2)  +  0(k3). 

Since  we  require  an  error  «  r,  we  set 

$k2(a3  +  a/Xz)  =  t 


giving 


k2 


6  r 

a(a2  +  l/X2)' 


Thus  w(r;X),  the  work  required  to  achieve  a  given  accuracy  t  using  Lax-WcndrolT  with  stepsize 
ratio  X,  is  given  by 


w(t;X) 


\W 
k 2 


= (Xa  + 1/Xa) 


a2W 

fir 
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Ji 


Wo  have  noL  yet  specified  X.  Choosing  X  to  minimize  w(r;  X)  gives  X  =  I  /a  and  the  (Minimum  work 
w(r)  is 


Mr)  = 


MW 

3  r 


for  unsplit  Lax-Wendroff. 


(2.8) 


Now  consider  the  split  method  (2.4a,b).  Let  Wa  be  the  work  required  to  apply  Lax- WendrofT  on 
the  slow  scale  and  Wl'f*v  the  work  required  to  computed  cxp(fc/l^<9I)f/JJl.  Set  WTSM  =  Wa  +  Wejxp. 
Typically  WTSM  »  W.  The  error  over  one  unit  of  time  for  the  split  scheme  is  bounded  by 


||((g/(A:/2)<3a(fc)Q/(/:/2))t/fc  -  oxp(/t<9x))  u|| 

<  ^||e8plit(fc)«  +  Ea{k)u  +  2/<;/(fc/2)u  +  0{k3) |) 


For  (2.4b),  Ef(k/ 2)  =  0.  From  (2.5a), 

||/?a(*HI  <  bKk2b3  +  h2b) 
=  i  k3(b3+b/\2). 


The  splitting  error  is  bounded  using  (2.2), 


!|£spiit(*H  <  ik3(Mb  +  ab*) 

«  kk3aH, 

although  it  may  he  much  smaller  for  some  problems.  Since  our  results  depend  very  much  on  the 
size  of  this  error,  we  will  suppose  for  now  that 

||0spiit(*HI  <  fk3o 


for  some  a,  so  that 

~\\ESf,i4k)u  +  E3(*M  <  \k\a  +  b3  +  6/X2). 
In  order  to  obtain  accuracy  r  we  must  take 


a  +  b3  +  6/X2 
so 

w(t)  X)  =  \WTSM/k2 

XA/TSM 

=  X(<r  +  63  +  (2-10) 

The  optimal  stepsize  ratio  X  now  depends  on  the  size  of  the  splitting  error  and  is  given  by 


X 


b 

a  +  b3 


so  that 


w(r)  =  \/6(ff  +  63) 


jyrsur 

3r 


for  the  time  split  method  (2.4a, b). 


(2.11) 
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(2.12) 


If  tr  <  63  (e.g.  when  A/  and  /t3  commute),  then  (2.11)  gives  X  as  1/4  ami 


w(t )  = 


b’WTSM 

St 


When  WTSM  »  W  this  is  better  than  (2.8)  by  a  factor  or  c2,  meaning  greatly  improved  efficiency. 
Note  that  when  a  =  0  the  only  error  incurred  is  the  error  in  using  Lax-  Wendrolf  on  the  slow  scale. 
From  our  previous  discussion  of  Lax- Wendrolf  it  is  dear  why  X  =  1/6  is  optimal  in  this  case. 

On  the  other  hand,  if  the  splitting  error  is  as  bad  as  (2.9)  indicates,  then  a  =  a2  6  and  X  rs  1/a 
in  (2.1 1)  giving 

,  s  abWT3M 

This  is  still  an  improvement  over  (2.8),  although  now  only  by  a  factor  of  e.  Note  that  now  X  is 
chosen  appropriate  to  the  fast  scale,  even  though  the  fast  part  of  the  problem  is  solved  exactly, 
in  order  to  reduce  the  error  due  to  splitting.  Indeed,  if  we  try  to  use  X  =  1/6  when  a  —  a2b,  we 
obtain  no  improvement  over  (2.8).  For  this  reason  it  is  advisable  to  always  use  small  time  steps 
with  the  time-split  method  (2,1a, b)  unless  Es,,iu (&)  <s  known  to  be  very  small,  in  which  case  even 
greater  efficiency  is  achieved  by  using  larger  limesteps. 

Now  consider  the  method  (2.1a,c)  where  Lax- Wendrolf  is  used  for  both  operators.  In  this  case 
WTSM  =  W„  +  mWf  where  W/  is  the  work  requried  to  apply  Lax- Wendrolf  on  the  fast  scale.  We 
are  assuming  that  Wj  W3  W .  We  will  take  m  «  Xa  as  suggested  in  (2.6).  Using  (2.5c)  we 
find  that 


rll^pli(.(£)«  +  E„{k)u  +  2Ef(k/2)u\\  <  $k2(r  +  63  +  6/X2  +  a3/m2  -I-  a/X2). 


We  then  obtain 


w(t;  X)  =  X(<r  +  63  6/X2  +  2 a/X2)  — a 

or 


X(a  +  63  +  2a/X2) W'  +-XaWr  for  (2,1a, c). 
or 


(2.13) 


The  optimal  X  now  depends  on  the  relation  between  Wj  and  W,  and  is  more  difficult  to  solve  for. 
We  will  discuss  three  possible  choices  of  X:  X  =  l /a,  X  =  1/6,  and  X  =  l/v/a6. 

When  X  =  1/a,  m  —  l  and  we  are  simply  alternating  between  LW(Af,  k)  and  LW{Aa,k).  In 
general  we  would  not  expect  this  to  be  any  more  efficient  than  using  the  unsplit  method  LW(A,  k). 
Indeed,  we  find  that 

W  -i-  W, 

w(t;  1/a)  «  ((<r  +  63)/a  +■  2a2) — - - 

6r 


2  W,  +  Wf 
1  St 


regardless  of  the  size  of  a.  This  is  better  than  (2.8)  only  if  Wa  +  Wj  <  W,  which  is  generally  not 
the  case  for  the  problems  we  arc  considering.  (Note  that  this  is  the  case,  however,  in  the  LOD 
method,  where  we  alternate  between  solving  onc-dimcnsiona!  implicit  problems  in  different  space 
dimensions.) 

For  X  =  1/6  increased  efficiency  is  possible  if  the  splitting  error  is  small.  From  (2.13), 


w{t ;  1/6)  =  (<r/6  +  262  +  2a6) 


W.  +  f  Wf 


6  r 


(<r/b  +  2a6) 


Wa  +  %W, 

6  T 
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Suppose  i.li at.  1  Vj  +  <  fV/  7  IF  for  some  7  <  1/2.  '1'lieu  IF,  +  £  Wy  —  ^  W  and  so 


-,iu 

/e(r;  1/6)  =  ( tia/lr  +  2a")  -  --. 


(2.11) 


This  is  belter  than  (2.8)  whenever  ci  <  alt2 /-f.  lti  this  ease  the  Lime-split  method  is  more  eflieient. 
For  example,  if  a  =  0,  (2.1  I)  is  belter  lhan  (2.8)  by  a  factor  of  7. 

Unfortunately,  if  a  =  u"6  as  in  (2.9),  then 

,  ,  ,L,  «*{w.  +  m) 

<"(t;  l/fa) » - 3^- - 

which  is  no  better  than  (2.8)  and  may  be  worse  if  Wj  >  rkF. 

Now  consider  an  intermediate  stepsize  ratio,  X  —  \j\fab.  From  (2.13), 

t  W3  +  yfiJbWf 

w(t;  l /'/alt)  - (o  +  2a" b) - 

6T 

~~  2r 

^  w,  +  \A/6  W7 

=  v^»  - * - 

regardless  of  the  size  of  tr.  This  is  better  than  (2.8)  if  Wj  -F  y/(W,  <  which  will  generally  be 
true. 

We  conclude  that  the  method  (2.4a,c)  is  more  efficient  when  X  is  chosen  correctly.  If  the 
splitting  error  is  known  to  be  small  then  X  =  t/6  can  be  used.  Otherwise  smaller  timesteps  should 
be  used,  e.g.  X  =  1  j\fah.  Very  small  timesteps,  X  =  1/a,  should  never  be  used. 

Here  we  have  not  dealt  with  the  advantages  of  the  split  scheme  resulting  from  the  possibility  of 
choosing  the  stepsize  ratio  on  each  scale  so  that  the  fc3  and  kh2  terms  in  each  of  the  Lax-  Wendroff 
errors  nearly  cancel  out.  When  this  can  be  done,  the  splitting  may  be  even  more  advantageous 
than  indicated  here. 

Block  triangular  systems.  Since  the  efficiency  of  the  split  scheme  is  limited  primarily  by 
the  splitting  error,  it  is  interesting  to  investigate  how  this  error  depends  on  the  coupling  between 
fast  and  slow  scales  in  a  simple  model  system.  Suppose  that  the  matrix  A  is  oT  the  Torm  (1-10) 
with  1 1 y\  1 2 1 1  ^  a  <  l  and  that  the  splitting  (1.1 1 )  is  used.  Here  A  |2  is  the  coupling  between  List 
and  slow  scales.  If  /li2  =  0,  the  problem  is  uncoupled  and  #Split (k)  =  0.  In  general,  from  (2.2), 

EspUi(k)=~  J  Ti^u(Mn^2-2/l,2/l22)  ^+0(Jki) 

Thus  ||Z?spijt(fc)u||  ak2 /2\c2.  The  efficiency  of  the  splitting  depends  on  the  size  of  a.  In  the 

notation  used  above,  we  have 

a  =  1/e,  6=1,  a  —  \aa2b. 

For  unsplit  Lax-WendrolT,  (2.8)  gives 

1  W 

W(r)=^-.  (2.15) 

The  time-split  method  (2.4a,b)  is  always  more  efficient  if  wc  choose 

X  «  (1  +  7Qa26)-1/2. 

For  example,  if  a  1  wc  should  use  X  m  2/a  =  2c  in  order  to  reduce  (2.15)  by  a  factor  of  c.  The 
maximum  efficiency  indicated  in  (2.12)  is  achievable  only  if  a  <  r2,  in  which  case  taking  X  =  1 
reduces  (2.15)  by  a  factor  or  c2. 


ifrMMrtifti 


Reducing  the  splitting  error.  For  block  triangular  systems  in  wliieli  A  |  j  is  not  sufficiently 
small,  it  is  possible  to  reduce  the;  coupling  through  a  change  of  variables  so  that  the  optimal 
elliciency  can  be  achieved.  A  change  of  variables  amounts  to  replacing  a  by  u  =  Ha  for  some 
nonsingular  matrix  II.  The  system  ut  =  Aux  then  becomes  ut  =  BAB~lux.  Clearly,  if  U  is 
chosen  to  be  the  eigenvector  matrix  of  A  then  the  problem  completely  decouples  into  independent 
scalar  equations.  We  are  seeking  something  less  expensive  which  only  decouples  the  fast  and  slow 
scales.  Thus  we  want  a  matrix  H  such  that 


BAB~l 


iClt  0  ' 

0  £22 


(2.16) 


with  ||Cu||  «  HC22II  «  1.  In  the  block  triangular  case,  it  suffices  to  consider  D  of  the  form 


B  = 

I 

B 12' 

,  B~l  = 

I 

B\2 

0 

/ 

0 

I 

Then 


BAB~l 


7 An  ~ \AuB\2  +  A\2  +  B^Aiz 

0  A22 


and  so  Z?i2  should  be  chosen  to  solve 


-Ay\B\2  —  B[  2A22  =  A12 
£ 


(2.17) 


in  order  to  completely  decouple  the  fast  and  slow  scales. 

In  the  present  context  solving  for  B t2  from  (2.17)  is  not  worthwhile.  In  order  to  achieve 
optimal  efficiency  we  need  only  reduce  the  coupling  by  one  or  two  factors  of  c.  Further  reductions 
do  not  gain  anything  once  the  Lax-Wendroff  errors  dominate.  This  suggests  taking 


B 12  =  cAti  A12 


(2.18) 


so  that 


where 


BAB-1  s 


A(,) 

/1|2 

A  22 


A\ 2  —  e/tn1/li2/l22. 


We  now  have  ||>1  iV ||  »  ca  provided  II  ^  l*  The  coupling  is  thus  reduced  by  a  factor  of 
e  through  the  use  of  a  very  simple  change  of  variables.  The  above  process  can  be  repeated  to 
obtain  additional  factors  of  e.  This  change  of  variables  has  been  suggested  by  Kreiss[9j  in  a  similar 
context. 

For  full  systems  of  the  form 


A  — 


A\2 

A22 


we  can  obtain  a  similar  reduction  in  the  size  of  both  off-diagonal  blocks  and  again  reduce  the 
splitting  error  by  several  orders  of  magnitude.  In  this  case  we  consider  B  of  the  form 


7  K 

7 

0 

7  +  KL 

K 

.0  /. 

L 

r 

L 

/ 
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II-  is  easy  l.o  verily  that  the  lower  corner  of  A  is  annihilal-eil  by  taking  L  to  satisfy 

—  LA  1 1  —  ■''I22  L  —  LA\<>L  +  A‘21  =0. 

£ 

Tile  matrix  K  can  then  be  chosen  as  before  to  remove  the  remaining  upper  corner.  This  results 
in  a  system  of  the  form  (2.1ft).  This  particular  transformation  is  discussed  more  completely  by 
O’Malley  and  Andorson[14j.  Again,  however,  we  are  not  interested  here  in  completely  annihilating 
the  corners,  but  rather  in  reducing  them  by  a  factor  off.  This  is  easily  accomplished  by  taking 

K  =  £ytr,%2 
L  =  -~cA‘2iAnl. 

Example  8.1  in  section  8  illustrates  the  use  of  the  change  of  variables  for  a  triangular  system. 

3.  Stability  of  the  time-split  method 

In  this  section  we  investigate  the  stability  of  the  time-split  method  when  applied  to  a  constant 
coellicient  problem  on  the  entire  real  line,  —  00  <  x  <  00  or,  alternatively,  on  a  finite  interval 
with  periodic  boundary  conditions.  When  Q2(k/ 2)  =  Qf(k),  as  is  true  for  the  splittings  (2.4),  for 
example,  Cauchy  stability  of  the  Strang  splitting  (1.8)  is  equivalent  to  stability  of  the  first  order 
splitting 

Un+l  =  Qa(k)Qf(k)Un.  (3.1) 

For  simplicity  we  restrict  our  attention  to  this  splitting,  and  set  Q(k)  =  Q,(k)Q  f(k). 

In  general  the  stability  of  Qa[k)  and  Qf(k)  does  not  imply  stability  of  Q(k).  Instead  stability 
must  be  checked  directly.  In  fact,  (3.1)  can  be  unstable  even  when  Q/(k)  and  Qa(k)  are  exact 
solution  operators  for  well-posed  hyperbolic  problems  as  the  following  example  shows. 

Example  3.1.  Let 

Af  = 

Then  the  problems  ut  =  A/ux,  ut  =  A„ux  and  ut  =  (Af  +  A3)ux  are  all  well-posed,  strictly 
hyperbolic  problems  for  any  value  of  the  parameter  p.  Let 

Qf(k)  —  exp(kAfdx),  Qa(k)  =  expJfcA,^*). 

and  let  G /(^,  fc/2)  and  Ga[£,  k)  be  the  corresponding  amplification  matrices.  For  the  exact  solution 
operators, 

GfU,k/ 2)  =  exp  (ik£Af) 

e*fcC  pi  sin  k£ 

0  e“*fc*  . 

Gsit.k)  =  exp(ik£Aa) 

cosfcf  tsinfc£ 
i sin  fc£  cos  k£ 

We  have  p(G/(<;,  fc/2))  =  p(Ga(£,k))  =  1  for  all  £  and  fc.  On  the  other  hand,  the  amplification 
matrix  G{£,  fc)  for  the  time-split  method  (3.1)  has  p(G(£,  fc))  =  1  for  all  £  and  fc  only  if  |^|  <  2. 
When  |/x|  >  2,  the  method  is  unstable.  Figure  3.1  shows  graphs  of  p((7(£,  fc))  for  p  =  5, 10. 

In  spite  of  this  example,  there  are  some  very  important  classes  of  splittings  for  which  the 
individual  stability  of  Q/(k)  and  Qa(k)  docs  imply  the  stability  of  Q(k).  It  is  useful  to  delineate 
such  classes,  since  the  stability  of  Qf{k)  and  Qa(k)  is  often  easy  to  determine,  whereas  the  stability 
of  Q(k)  may  be  quite  tedious  to  determine  directly. 
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Figure  3.1.  Spectral  radius  of  the  amplification  matrix  C(£,  k )  of  example  3.1  for  p  =  5, 10,  as 
a  function  of  £k  between  —x  and  x. 


Block  triangular  systems.  One  such  case  is  the  block  triangular  system  of  equations 


u 

Mil 

A12 

XL 

V 

t 

0 

A'i  2 

V 

with  the  splitting  (1.11).  The  solution  v  docs  not  depend  on  u.  In  solving  for  u,  the  computed  v 
enters  essentially  as  a  forcing  function.  The  schemes  Qa{k)  and  Qf{k)  will  be  of  the  form 


<?.(*)  = 


I  Qu{k) 
[0  Qinik) 


QfW  = 


Qu(k)  0 

0  I 


(3.2) 


Suppose  that  Qi  i(fc)  and  Q22(k)  arc  stable  schemes,  and  in  particular,  that  there  exist  a  norm  ||  •  || 
and  a  constant  a  >  0  such  that 


||Qu(fc)||  <  1  +  ak  for  all  k  sufficiently  small.  (3.3) 

AH  of  the  following  estimates  will  be  in  this  norm.  We  also  suppose  that 

||Q12(*)K||  <  kM\\D+V\\  (3.1) 

for  some  constant  M.  These  assumptions  arc  satisfied  for  the  methods  (2.1)  provided  the  Lax- 
WondrolT  operators  arc  stable.  Wc  then  have  the  following  theorem. 


13 


THEOREM  .‘{.I.  Suppose  Qf(k )  .iiul  (J„ [k)  are  stable  schemes  as  above.  Then  the  split,  scheme 
(Ji(k)O f(k)  is  stable  for  smooth  initial  data  l/fl.  More  precisely,  wo  obtain  bounds  for  the  solution 
which  depend  on  a  discrete  Sobolev  norm  of  the  initial  data, 

IKm  <  /vt(||//°1|  +  ||/>+v'0||)  (:*.5a) 

\\Vn\\  <  /<t||^0||  (:*-r>b) 


(3.6a) 

(3.6b) 


for  nk  <  T.  Here  K?  and  Kt  are  constants  depending  only  on  the  fixed  lime  T. 

Proof.  When  the  lull  scheme  (/n  +  1  =  Q,(k)Qf(k)Un  is  written  out  we  obtain 

Un+'  =  Qu{k)Un  +  Qu(k)Qv2(k)Vn 
K"+'  =  Q22{k)Vn 

The  bound  (3.5b)  follows  immediately  from  (3.6h)  and  the  stability  of  Q22(k).  Moreover,  by 
linearity,  an  identical  bound  holds  for  the  linear  combination  of  solutions  D+Vn,  i.e., 

\\D+Vn\\  <  kT\\D+V°\\. 

Using  this  together  with  (3.1)  in  (3.6a)  gives 

\\Un+'\\  <  IK^i i(^)ll(ll^/n|l  +  kM Kt\\D+V°\\). 

When  iterated  n  times  this  gives 

l|t/n||  <  HQii(*)iril^°ll  +  kMI<T{\\Qn(k)\\n-1  +  HQiiWII"-2 

+  ---  +  ||Q.i(fc)||  +  l)||D+K°||. 

IJy  (3.3),  ||gu(A:)ir  <  (l  +  afc)n  <  eaT  if  nk  <  T.  Using  this  in  (3.7)  gives 

\\un\\  <  fiaT(||g°||  +  TMKt\\D+Vq\\) 
for  nk  <  T,  which  is  of  the  desired  form  (3.5a).  | 


(3.7) 


Simultaneously  normalizable  matrices.  Stability  also  follows  directly  when  Af  and  A, 
arc  normal  matrices  (a  normal  matrix  is  one  which  commutes  with  its  transpose).  This  includes, 
for  example,  symmetric  matrices  and  scalar  problems.  In  Tact,  it  suffices  that  Af  and  Aa  be 
simultaneously  normalizable,  i.e.,  that  there  exist  some  nonsingular  matrix  S  such  that  S/ts5~l 
and  SAfS~l  are  both  normal.  Thus,  the  case  of  simultaneously  diagonalizablc  Af  and  Aa  is  also 
covered.  This  is  a  consequence  of  the  following,  even  more  general,  theorem. 

THEOREM  3.2.  Let  A\,A2,  ...,Am  be  constant  matrices.  Approximate  each  solution  operator 
cxp(kjAjdx)  by  some  operator  Qj[kj)  with  amplification  matrix  Gy(£).  Suppose  there  exists  some 
norm  ||  •  ||  for  which 

l|C,-(flll<l  Vf,  j  =  l,2,...,m.  (3.3) 

Then  the  scheme 

Un+X  =Qi(kl)(h(k2).--Qm(km)Un  (3.1) 

is  stable. 


Proof  Let  G(£)  =  G’i^K'sU)’  ”Gm{€)-  Then  powers  of  (7(£)  are  uniformly  bounded  in  the 
norm  ||  ■  ||  since 


\\cnm  <  nc(oir 

<  (iici(oii-  •  •iK-m(oii)n 

<  i. 


It  follows  that  (3.1)  is  stable.  | 


(.'OUOI.I.AUY  Suppose  liter e  exists  some  non.smgu/ar  matrix  S  such  that  SAjS'  1  is  normal 
tor  j  —  1,2,  . ..,  m  ami  that  Liu:  ami>Hlicaiiou  mat  rices  (l}(£)  satisfy 

0  P((:A0)  <  1  v£>  J  =  l,2,...,m 

ii)  SCI j(£).9— 1  is  also  normal  for  all  £,  j  =  l,2,...,m 

T/ien  t/ie  scheme  (3.1)  is  stable. 

Remark:  Condition  (3.5tt)  is  satisfied  if  Qj(k})  is  the  exact  solution  operator  or  one  or  more 
steps  of  Lax-Wendroff. 

Proof.  Since  the  2-norrn  of  a  normal  matrix  is  equal  to  its  spectral  radius,  conditions  (3.5) 

give 

\\sc,{t)s-%  =  =  P(CAO)  <  I- 

It  follows  that  the  hypothesis  of  Theorem  3.2  is  satisfied  in  the  norm  ||  •  ||  defined  by 

||/1||  =  ||S/1S-1||2. 


This  completes  the  proof.  | 


4.  The  Leapfrog  Duhamel  method. 

As  mentioned  in  the  introduction,  the  time-split  method  does  not  immediately  lend  itself  to  use 
with  multi-level  difference  schemes.  We  now  present  a  new  method  with  the  same  basic  philosophy 
as  the  time-split  method  but  which  uses  leapfrog  on  the  slow  time  scale. 

Using  Duhamcl's  principle  (i.e.,  variation  of  parameters)  we  can  write  the  solution  to  (1.3)  as 

rt+k 

u(x,  t  +  k)  =  exp(2fc/t j rJx)u(x ,  t  —  k)  +  /  exp((t  +  k  —  r)A  fDx)A3ux(x,  r)  dr. 

Jt-k 

If  we  now  approximate  the  integral  by  the  midpoint  rule  wc  obtain 

u(x,  t  +  k)  «  exp(2kA/dx)u(x,t  —  k)  +  2kexp(kAfdx)A,ux(x,t) 

=  exp(fc/t /dr)[cxp(fc/t  fdx)u(x,  t  —  k)  +  2kA„ux(x,  <)]. 

Replacing  ux(x,t)  by  the  standard  centered  difference  operator  and  approximating  exp(kA/Dx)  by 
Q/(k)  gives  the  Leapfrog  Duhamel  method, 

tc+l  =  Q/MQ/WZ -l  +  ju.(tc+l  -  fC_,)].  (4.1) 

The  term  inside  the  brackets  is  essentially  leapfrog  for  the  problem  ut  =  A,ux  since  Q/(k)Un~l  ss 
cxp(  —  kA,t)x)Un.  If  Qf(k)  is  an  0(k3)  approximation  to  exp(kAfOx)  then  ('1.1)  provides  an  0(k3) 
accurate  approximate  solution,  even  for  noncommuting  Ay  and  Aa.  This  will  be  shown  in  section 
5  where  the  method  is  analyzed  in  more  detail. 

Wc  pay  a  price  for  using  a  scheme  involving  three  time  levels,  since  (LI)  requires  two  applica¬ 
tions  of  the  operator  Q/(k).  One  of  these  is  needed  only  to  provide  the  proper  values  at  time  n—  1. 
Nevertheless,  this  method  may  be  useful,  particularly  in  cases  where  cxp(ltyl/dI)  is  known  exactly 
and  thus  is  easy  to  apply. 
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5.  Accuracy  of  the  Leapfrog  Duhamel  method 

The  Leapfrog  Duhamel  scheme  can  he  analyzed  hi  terms  of  the  error  in  the  midpoint  rule 
directly  from  its  derivation.  We  prefer  to  build  upon  the  results  in  section  2  by  rewriting  Leapfrog 
Duhamel  as  a  splitting. 

First  consider  the  standard  leapfrog  scheme  on  Ut  =  Aaux, 

Un+i  =  Un~l  +  2  kAaD0Un. 

The  truncation  error  is  given  by 

[u(x,  t  —  k)  +  2kAaD0u(x,  t)j  —  u(x,  t  +  fc) 

=  [{[  +  2 kA,Do  exp (kAadx))  —  exp(2fc/lsdI)]u(i,  t  —  k). 

For  conformity  with  section  2,  we  define  the  operator  Qa(2k)  by 

Qa(2k)  =  [  +  2kAaD0oxp{kAadx). 

Note  that  this  is  not  the  actual  finite  difference  operator  for  leapfrog,  since  in  general  Un  is  not 
exactly  equal  to  exp(fcyl3dI)t/n_l ,  but  it  is  the  proper  operator  for  computing  the  local  truncation 
error,  in  which  Un~ 1  and  Un  are  replaced  by  the  true  solution  values.  We  can  now  define  the 
truncation  error  operator  for  leapfrog  on  stepsize  k  by 

l'^F(k)  =  Qa( 2k)  -  exp(2 kA.dx)  =  0(k3).  (5.1) 

The  Leapfrog  Duhamel  scheme  is 

Un+l  =  Qf{k)(Qf(k)Un-'  +  2kAaD0Un)  (5.2) 

whore  Qf(k)  is  some  approximation  to  cxp(kAjdx)  with  error  operator 

Bf{k)  =  Qj(k)  -  cxp(kAfdx)  =  0(k3). 

To  obtain  the  truncation  error  for  Leapfrog  Duhamel  we  replace  Un~l  and  Un  by  u(x,t  —  k)  and 
u(x,  t)  in  (5.2).  The  right  hand  side  then  becomes 

Q/{k)[l  +  2kAaDo  exp(fc/td*)Qj '(fc)]Q/(fc)tt(z,  t  —  k) 

=  Q/(k)[f  +  2kAaD0  cxp(fc/ts5*)  +  2kAaD0(exp(kAdz)QJl(k) 

-  cxp{kAadx))]Q  f{k)u{x,  t  -  k)  (5.3) 

=  [Qf(k)Qa(2k)Qf{k)  +  2kQf{k)AaD0(cxp(kAdx) 

—  exp (kAadx)  cxp(A:/l/9I)  —  cxp(kAdx)Ef(k))]u(x,  t  —  k). 

Thus  the  Leapfrog  Duhamel  operator  can  be  viewed  as  a  splitting  of  the  form  Q f{k)Q a{2k)Q f(k) 
plus  some  additional  error  terms  which  arc  0(k3).  Let  Firap\it[k)  denote  the  error  operator  for  the 
first  order  accurate  splitting 

KfPlii(*)  =  exp(Mt>*)  -  cxp(kAadx)  cxp(kAfdx)  =  0(k2). 


Observing  Dial 


exp(kAai)x)lif(k)  =  <){k'), 

Qf(k)  =  /  +  0{k), 

Do  =  <9*  +  D(k“), 


the  operator  in  (5.3)  becomes 

Qf[k)Q.{2k)Qf{k)  +  2  kAadxEiplit(k)  +  ()(kA). 

Using  (2.3)  we  obtain  an  expression  for  the  truncation  error  operator  for  Leapfrog  Duhamel, 

liLFD(k)  =  (Qf(k)Qa(2k)Qf(k)  +  2  kAadxKlpVll(k)  +  0(k 4)) 

—  cxp(‘2kAdx) 

=  E.Piti(2*)  +  U^F(k)  +  2  lif{k)  +  2kA,<)J<;'lit,iit(k)  +  0{kA). 


For  Aa  and  Af  constant  we  have 


Elp iu(*)  =  \k2(AfA.  -  AaAf)dl  +  0(k3) 

so 

£.Piit(2*)  +  2kAadxElpni(k) 

=  -  k k\A)Aa  -  2AfAaAf  +  AaA) 

+  A]Af  +  AaAfAa  -  2 AfAi)d\. 

The  splitting  error  in  Leapfrog  Duhamel  is  thus  roughly  8  times  as  large  as  the  corresponding  error 
in  the  time  split  method  with  Lax-Wendroff.  The  work  comparisons  of  section  2  can  be  repeated 
for  Leapfrog  Duhamel  with  similar  results. 


6.  Stability  of  the  Leapfrog  Duhamel  method 

At  present  the  stability  analysis  for  Leapfrog  Duhamel  covers  only  the  case  in  which  Af  and 
Aa  are  simultaneously  diagonalizable, 

XA/X~l  =  Mf,  XAaX~l  =  M, 

where  Mf  and  Ma  are  diagonalizable  matrices.  We  assume  that  Qf(k)  is  stable  and  is  also 
diagonalized  by  X.  This  is  true  for  Qf{k)  —  exp (kAfdx)  or  for  Qf(k)  =  ( LW(Af,k/m))Tn  with 
p(Af)k/mh  <  t.  Let  qf(k)  be  a  single  diagonal  clement  of  XQf(k)X~l  and  pa  a  diagonal  element 
of  Ma.  It  suffices  to  consider  the  scalar  equation 


Un+l  =  q){k)Un~X  +  2  kqf(k)paD0Un.  (6.1) 

Let  <7/(0  be  the  amplification  factor  corresponding  to  qf(k).  By  assumption,  |<7/(0I  <  1  f°r  £• 

THEOREM  6.1.  Suppose  ]X/xa{  <  l,  where  X  =  k/h.  Then  the  amplification  factor  </(£)  for 
the  scheme  (6.1)  satisfies 

IjKOI  =  lff/(0l- 
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t'rool'.  Tin'  amplification  factor  is  derived  by  letting 


=  <jntiwimh 

in  (6.1).  Wo  obtain  the  equation 

</(£)  =  ff/(C)ff-,(0  +  ii\gf(S)n»s\nth 

which  can  be  rewritten  as 

(3(0fl/‘(0)2  -2iX/i„sin£/ii ■(</(£)<// ‘(0)“  1  =  0. 

Solving  this  quadratic  equation  yields 

»U)0/l(f)  =  s,n  ± 

If  |X/i3|  <  l,  the  square  root  is  real  and  so 

I  <jti)9jl(Z)\2  =  l 

and  hence 

m\  =  iff/(oi 

as  claimed.  | 

Note  that  when  the  exact  solution  operator  is  used  for  q/(k)  we  have  |j//(OI  —  1  and  hence 
lfl(OI  —  1  f°r  nil  £.  In  this  ease  Leapfrog  Duhamci  is  nondissipativc. 

7.  Boundary  data  for  the  intermediate  solutions. 

For  general  initial  boundary  value  problems  we  must  be  able  to. generate  the  appropriate 
boundary  values  for  the  intermediate  solutions  which  arise  in  the  use  of  a  split  scheme.  We  have 
developed  a  general  methodology  for  defining  the  proper  boundary  data  which  will  be  illustrated 
here  for  constant  coefficient  problems  at  an  inflow  boundary.  More  general  problems  can  also  be 
handled,  as  will  be  reported  on  elsewhere.  The  procedure  will  be  demonstrated  for  the  time-split 
method  (2.4a,b),  but  can  also  be  used  for  the  other  methods  previously  described. 

■First  consider  the  scalar  problem 

=  ~(1  +  f )«,  x  >  0,  t  >  0 

«(0,  t)  =  g(t)  t>  0  [  ] 

with  the  splitting 

A{  =  -1,  A,  =  -€. 

Take  k  =  2h  and  use  the  method  of  characteristics  solution  for  Aj  and  Lax-Wcndroff  on  .1,  . 
There  is  no  need  to  use  a  Strang-type  splitting,  since  the  operators  commute,  and  thus  the  split 
scheme  is  simply 

=  m  =  2,3,... 

UVl  =  V'm~  -  lC-i)  +  2c2(f/^+,  -  2U'm  +  Ul_x)  (7-2) 

m  =  1,2,... 
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Tl»'  value  of  f./JJ  + 1  is  given  by  I  la-  boundary  conditions, 

=<j(tn+l). 

For  the  splitting  (7.2)  wc  must  also  provide  f/,*(  and  IJ\.  In  general  with  k  =  ph  for  some  integer 
p  >  1,  we  would  need  to  supply  IJ\, . . .,  //*_,. 

In  order  to  generate  boundary  data  wc  consider  (J’m  as  an  approximation  to  ii‘(im,(n+i) 
where  the  continuous  function  u’[x,t)  satisfies 


ut  —  ui  x  —  0|  t  >  tn  (7  3) 

u'(x,  t„)  =  u(x,  tn)  x  >  0. 

Then,  using  the  differential  equations  governin'.'  a  and  u’ ,  we  can  express  U'Q  and  U\  in  terms  of 
g(t).  Consider  U0.  We  want 

(/'  =  n*(0,S  *  *  ■ 

=  «*(0,  t,x}-  *  -• .&,  !„)  +  $k2u‘u[ 0,  «„)  +  •••  (7.4) 

=  uH),  t  .  \  •  tn)  +  |l:2u’t(0,  tn)  +  •  •  • 


Here  we  used  17.3)  to  express  in  terms  of  a*.  But  since  u'(x , <„)  =  u[x,tn)  for  all  x,  this 
relation  can  be  differentiated  with  respect  to  x,  giving  u*(x,  <„)  =  uz(x,tn)  and  similarly  for  higher 
derivatives.  So  (7.4)  becomes 

U o  =  u{ 0,  <„)  -  kux( 0,  tn)  +  \k2uxx{ 0,  /„)  +  •••. 


We  can  now  use  the  original  equations  (7.1)  governing  u  to  rewrite  this  in  terms  of  ^derivatives  of 
u.  Since 


j  >  0 


wc  obtain 

f/0  =  w(0,  tn)  +  wt(0,  tn)  +  j^7)2a«(0,  <n)  +  •  •  • 

=  u(0,  tn  +  A:/(l  +  <0)  (7-5) 

=  9(tn  +  k/{l  +  e)). 

This  is  the  desired  boundary  data. 

For  such  a  simple  example  it  is  easy  to  verify  that  this  is  the  correct  boundary  value.  According 
to  the  scheme  (7.2)  we  would  really  like 


Uo  =  Ul2  =  u(-2h,tn). 


Of  course  u  is  not  officially  defined  for  x  <  0,  but  using  the  differential  equation  (7.1)  it  can  easily 
be  extended  from  the  boundary.  Since  (7.1)  has  characteristics  with  slope  1/(1  +  e),  wc  find  that 


u(—2h,  tn)  =  u(0,  tn  +  2hj[\  +  e))  =  g(tn  +  */(!  +  t)) 

exactly  as  in  (7.5)  . 

Wc  can  compute  U ,  in  the  same  manner.  Wc  want 

U\  =  u*(M»+t) 

=  u*(0,  tn+ 1/2)  where  tn+i/a  =  tn  +  k/ 2. 
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Wo  now  procwd  as  before, 


U  |  —  it  (l),  /„)  +  .j  A  ut  (0,  tn)  +  ljk~utt( 0,  <„)  +  ••• 

=  u(0,t„)-  \kux(Q,tn)+  J-fc2«XI(0,  tn)  +  ■■■ 

=  «(0,  tn)  +  ^  ut(0,<n)  +  £{  T^)«tt(0,  <n)  +  *  ■  • 
=  </(<«  +  V  A/(l  +  r)). 


To  summarize  our  procedure,  we  switched  from  ^-derivatives  of  u  to  z-dorivatives  of  u  . 
Since  those  were  evaluated  at  time  tn,  they  were  identical  to  the  corresponding  z-dorivatives  of 
u.  VVe  then  switched  back  to  ^-derivatives  of  u  along  the  boundary,  which  allowed  us  to  use  the 
known  boundary  conditions  for  u.  Clearly  this  procedure  will  not  work  so  neatly  when  we  deal 
with  variable  coefficients,  systems  of  equations,  or  inflow-outflow  boundaries.  Nonetheless,  these 
same  ideas,  combined  with  a  little  ingenuity,  lead  to  sufficiently  accurate  approximate  boundary 
conditions  for  a  wide  variety  of  problems. 

Constant  coefficient  systems.  Next  consider  the  system  of  equations 


ut  =  Aux  =  (A f  +  Aa)ux  x  >  0,  t  >  0 
«(0,  t)  =  g{t)  t  >  0. 


We  assume  that  A  and  Af  have  strictly  negative  eigenvalues,  fn  general  A/  and  Aa  do  not  commute, 
so  we  will  have  to  use  a  Strang-type  splitting.  There  will  be  at  least  two  intermediate  solutions, 
say 

U‘  as  cxp(±kAfdx)Un 

(/**  «  exp(kAaDx)cxp(}}kAfdx)Un. 

Of  course  there  may  be  many  more  if  cxp(^kAfdx)  is  itscIT  approximated  by  several  steps  of  Lax- 
WendrolT,  but  they  can  be  handled  similarly.  The  general  principle  should  be  clear  from  considering 
(7’8)- 

Again  let  u  (x,t)  be  a  continuous  function  satisfying 


u*t  =  Afu‘x  x  >  0,  t  >  tn 

u‘(x,tn)  =  u(x,tn)  x  >  0. 


(7.9) 


We  then  want 

Uq  =  u  (0,  1/2) 

=  u*(0,tn)  +  %kut[ 0,  «n)  +  jjk2U*tt(0, *„)  +  ••• 

=  w(0  ,t„)  +  £fc/i/ux(0,  tn)  +  gfc2/t2uxx(0,  tn)  +  •••  (7.10) 

=  u(0,tn)  +  £fc/l//l_1Ut(0,fn)  +  A&2/t^.A_2Utt(0,  *„)  +  ••• 

=  9{tn)  +  \kAfA~x  g'(tn)  +  %k2A2A~2g"(tn)  +  • 

We  assume  that  the  boundary  is  non-characteristic  so  that  A  is  invertible.  In  general  U$  must 
now  be  approximated  by  the  first  few  terms  of  (7.10).  If  we  keep  only  the  first  two  terms  wc  will 
have  boundary  data  with  0(k2)  errors.  This  is  sufficient  to  retain  the  0(k2)  global  accuracy  of 
Lax-Wendroff  (sec  Gustafsson(7)).  It  may,  however,  increase  the  error  constant  considerably  and 
partly  offset  the  benefit  obtained  by  using  the  split  scheme.  Consider,  for  example,  a  case  in  which 
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||.ty||  1  nml  ||.l.,||  ~  (.  Ill  this  ease  the  error  Hr*  M  (fc)u  is  like  <  k6  al.  most,  and  Hie  resulting 

global  error,  assuming  it  is  smooth,  will  be  like  <k~.  In  order  to  achieve  the  same  aceiiraey  in  the 
boundary  data  we  will  have  to  include  the  third  term  of  (7.10)  as  well  (or  at  least  its  dominant 
part).  In  some  such  cases  it  happens  that 

AjfA~’  =  I  4-  (){c)  for  j  =  1 , 2, . . .. 

We  can  then  retain  0(tk2)  accuracy  simply  by  taking 

Uo  ~  (j{tn+  i/a)  +  1  —  7)t/'(tn). 

This  will  be  illustrated  in  Example  7.1  . 

Now  to  find  boundary  values  to  U*  .  The  easiest  way  to  proceed  is  to  note  that 

U’*  =  cxp(—  kkAfOx)Un+l 

which  prompts  us  to  define  u  '(x,t)  ns  the  continuous  solution  to 

(®» 0  =  (*»  0  *  ^  ^  ^  ^n+l 

u"(x,  <n+1)  =  u(x,  t„+i)  x  >  0. 


We  now  solve  this  backwards  in  time  for 


U0  — tt  (0,  <n+t/2). 


Proceeding  as  in  the  derivation  of  (7.10)  wc  obtain 


Uo  =  a{tn+\)  ~  r^kAfA  V(<n  +  |)+  zk*A)A  V'(*»+l)  + 

ffUn  +  1/2)  -  ^(AfA~l  -  f)(j’[tn+  l)' 


Example  7.1  Consider 


«1  f— 1  £i  11V 

=  „  0  <  x  <  1,  t  >  0 

U2  -2JM, 

u(x,  0)  =  f(x)  0  <  x  <  1 
u(0,  f)  =  g(t)  t  >  0 

where  u  =  (u,  u)T.  We  have  chosen  a  strip  problem  to  illustrate  that  outflow  boundaries  are 
frequently  trivial  to  handle  with  a  split  method.  Take 

,  r-i  01  _ro  t , 

0  — 2J  Us  0 

For  this  problem  the  splitting  error  is 


£»Plit(fc)  —  g*3  l£a  ^  dl • 


If  we  list;  the  time-split  method  (2.1:i,l>)  then,  according  to  (2.1 1 ),  the  optimal  steps!  ze  ratio  is 


where  c  =  max|ty|.  For  k  =  2h  and  h  =  1/M,  (2.4a,b)  becomes: 

m  =  1,2,...,  M 
Vl  =  V^-2,  m  =  2,3,...,M 


Um  =  LW(Aa,k)Um,  rn  ■=  1, 2, . . M  —  1 


f/n+‘ 


f/ 


0  =  £/(<n-«-i) 

n-f !  _ 

m—  1  > 


=  u 


m  =  1,2 M 


^m+'  =  C-2>  m  =  2,3,...,M 


Notice  that  no  boundary  conditions  whatsoever  need  to  be  specified  at  the  outflow  boundary  x  =  1. 
On  the  inflow  side  we  still  need  to  specify  U0,  V\,  IJQ  ,  and  V,  .  For  ths  problem, 

Ma-1  =  1  +  3" 

7  (2-fi€2)2L  12c2  4  +  4eit2 

=  l  +  0(c). 

and  we  can  retain  0(ek 2)  accuracy  by  taking 


=  ff(*«+i/«)  +  hk(A,A~y  -  l)g'(tn) 
k 


=  3(tn+l/2)  + 


f |t2  Cl 
2t2  c i c2 


»'(««) 


(7.14) 


2(2  —  c  i  c  2 )  1 

Similarly  we  use 

VO  =  g(*n+1/2)  -  -  ')»'(*«  +  .)• 

In  order  to  implement  the  split  scheme,  we  also  need  V\  and  K7+l-  We  want  V\  = 
v*(/i,<n-t-i/2)  —  v* (0,  tn+\/\)  and  so  the  appropriate  value  comes  from  the  second  equation  of 

u*{0,tn+l/4)  as  g{tn+ 1/4)  +  \k(AfA~l  - 


i.e., 

V\  =  <72(«„+l/4)  +  4(r-«.«,)(2ggg>(f")  +  f‘  €^i(tn)) 
where  <7  =  (j;i,<72)r.  Similarly, 

Vr?+1  =  ff8(tn+3/4)  “  4(2jt|ta)(2<:2g>l(^n4-l)  +  C 1  C2ff2(tn+ 1 )). 

Computations  confirm  that  these  boundary  conditions  preserve  0(tk2)  global  accuracy  in  the 
split  scheme.  Actually,  for  this  particular  example  with  k  =  2 h,  even  greater  accuracy  can  be 
achieved.  Computing  Fi,(k)  from  (2.5a)  shows  that  the  0(ek 3)  terms  exactly  cancel  the  0(ek 
terms  in  £^apiit(^)  t  and  that  the  total  truncation  error  I£TSM(k)u  is  actually  0(e2fc3),  giving  0(e2k2) 
global  accuracy.  Higher  order  boundary  condit’ons  can  be  derived  which  maintain  this  accuracy, 
but  this  cancellation  of  errors  is  a  fluke  which  docs  not  occur  in  general. 
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Stability  for  the  initial  boundary  value  problem.  The  boundary  approximations  derived 
here  all  depend  only  on  the  given  boundary  function  </(t)  and  its  derivatives.  Suppose  the  time- 
split  method  used  in  the  interior  is  Cauchy  stable.  Then  the  stability  of  the  result i u p;  scheme 
for  the  initial  boundary  value  problem  follows  directly  from  the  theory  of  Ciisfafsson,  Kreiss  and 
Sundstrbm[8],  if  we  modify  their  stability  deliniLion  3.3  by  using  an  appropriate  Sobolev  norm  of 
the  boundary  data  on  the  right-hand  side. 


8.  Computational  results. 

In  this  section  we  give  various  examples  of  splittings  and  present  the  results  of  some  numerical 
experiments.  The  first  example  is  a  2  X  2  upper  triangular  system  of  the  form  (1.11).  We 
demonstrate  the  effects  of  the  splitting  error  and  its  reduction  by  the  use  of  a  simple  change 
of  variables  as  discussed  in  section  2.  1 

The  second  example  is  a  variable  coefiic'n  nt  scalar  equation  in  which  the  coefficient  lias  small 
variations  around  some  large  mean  value.  We  give  an  expression  for  the  splitting  error  in  such 
pro  blems. 

In  example  8.3  we  consider  the  one-dimensional  shallow  water  equations.  In  some  cases  this 
system  can  bo  broken  up  into  a  constant  fast  part  and  a  quasilinear  slow  part  in  conservation  form. 

Example  8.1.  This  problem  is  designed  to  illustrate  the  effects  of  the  splitting  error.  Consider 


ut 


to 

0 


for  0  <  x  <  1,  t  >  0 


(8.1) 


with  initial  conditions 

«i(z,0)  =  0)  =  e-,00(I-,/2)a 

and  periodic  boundary  conditions 

u3(0,  t)  =  «,(M)  <>0,  j=l,2. 

Figure  8.1  shows  the  results  after  236  time  steps  using  Lax-WendrolT  with  h  =  1/50  and  k  =  h/  10 
on  the  unsplit  problem.  Figure  8.2  shows  the  results  based  on  the  splitting 


Af 


10  0 
0  0 


A. 


0  f 
0  1 


We  used  k  —  h  =  1/50  with 

Q.(k)  =  LW(A„  k),  Qf(kf  2)  =  {LW(Af,  k/ 10))5. 


In  this  case  E3(k)  =  Ef(k/ 2)  =  0  by  a  judicious  choice  of  k/h  and  m.  The  second  component  U2 
is  computed  exactly  and  the  errors  in  u\  arc  due  entirely  to  the  splitting  error. 

If  the  change  of  variables  suggested  in  (2.18)  is  applied  twice  to  (8.1)  with  f  =  0.1,  wc  obtain 
the  new  variable 

u,  =  Ui  —  («  +  f2)«2  =  U\  —  0.11U2 


and  (8.1)  becomes 


«2 


10  0.0! 

.0  1 

«2 

it 
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II  we  solve  I  his  system  with  the  same  split  scheme  as  before  and  then  t  raiisronn  link  to  (.lie  original 
variables  by  u(  =  a )  +  O.llu-j,  the  errors  in  ii\  are  reduced  to  0(1(1  *)  as  mt:i  is  figure  8. lb 

The  Leapfrog  Dohamel  scheme  can  be  applied  to  this  system  with  similar  results.  'I’lie  saimt 
change  of  variables  ran  clearly  be  used  to  reduce  the  splitting  error  in  this  scheme  as  well. 

Example  8.2.  For  problems  of  the  form 

ut  =  (a  +  n(r.))ux 

with  a  constant  and  |«(x)|  <£  |a|,  the  splitting  error  operator  corresponding  to  A/  =  a,  A,  =  u(x) 
is 

Eapt,i{k)  =  cxp(  ^kaOz)  exp(Am(z)<Tx)  cxp(  1  kat)x)  -  cxp(fc(a  +  a(x))i>T) 

=  -  («'(x))2)/)x  +  0(fc<). 

For  the  Leapfrog  Duhamcl  scheme  the  splitting  error  is 

li*M  2k)  +  2ka(x.)t)xliiplil(k) 

—  ^spiit(2^)  +  2  ka(x)()z(^k'a(\i'(x)<)x  +  0(&3)) 

=  —  ^A:3a[(2a  +  a(z))«"(z)  -  'l(t*'(z))2  —  .‘la(x)r*'(z)<i>I]dI  +  OfA:3). 

The  Lax-Wendroff  and  leapfrog  errors  on  ut  =  a(x)ux  are  respectively 

^W(k)  =  -^M*)[«2(z)r?3  +  3rr(z)a'(z)rT2  +  ((o'(z))2  +  a(z)«"(xpz] 

+  \bh2n(x)<)l  +  0(k'{) 

0 

and 

EL/{k)  =  2  E^w(k)  +  ()(kA). 

For  the  lest  problem 

ut  =  (l  +  0.1  sin(27rz))«x  on  [0, 1] 
u(x,  0)  =  sin(l7rz)  0  <  z  <  I 
u(0,t)  = 

a  comparison  of  the  errors  shows  that  the  splitting  error  for  either  scheme  with  k  =  \h  should  be 
of  roughly  the  same  size  as  E,(k)  and  considerably  sinallc-  than  the  error  for  the  unsplit  operator 
with  the  same  spatial  step  and  reduced  time  step  k  =  hj 2.  Thus  vve  expect  the  split  scheme  with 
the  true  solution  operator  used  on  ut  =  ux  to  be  more  accurate  than  the  unsplit  scheme.  This  is 
confirmed  by  the  computational  results  in  Table  8.1.  Note  that  in  this  case  the  improver!  accuracy 
was  obtained  using  only  about  one  eighth  the  work  required  for  the  unsplil  scheme. 

If  Lax- WendrolT  is  used  on  the  fast  scale,  Qf(k/ 2)  —  [LW{Af,  k/ 8))3,  the  corresponding  error 
2Ef{k/2)  is  roughly  the  same  size  as  the  error  in  the  unsplit  scheme.  This  error  dominates  in  the 
resulting  split  scheme  and  so  we  get  roughly  the  same  accuracy  as  in  the  unsplit  scheme.  This  is 
also  illustrated  in  Table  8.1. 

Example  8.3.  The  onc-dimcnsionai  shallow  water  equations  can  be  written  as 


V 

V 

l 

V 

w 

t 

[<t> 

V 

0 
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where  n(x ,  t)  is  l.  In'  velocity  and  <j>  =  tjh  with  h(x,t)  I  In:  height  and  <j  flu:  gravitational  constant. 
Typically  /)  =  i,i  +  t)  where  0  is  constant  and 

|0'(M)|  «  |t/>| 

\  »(x>  01  <  !<£!• 


With  the  change  of  variables 

u{x,  t)  =  <j)  v(x,  t) 

the  system  (8.2)  becomes 

u 

4> 

The  natural  splitting  is  then 

0  l 
I  0. 


Aa{u,<f>)  =  -tf> 


u  ,) 

<t>'  u 


=  -4> 


-1/2 


.  “  4> 

*  +  ui 


We  have  ||/13||  ||/1/||.  The  matrix  Af  is  constant  and  the  method  of  characteristics  can  easily 

be  used  for  Qj(k/‘2)  .  Furthermore,  the  problem  on  the  slow  scale  can  be  written  in  conservation 
form.  Since  4>x  —  <j>'x  we  have 


/ta 


;l  -  (-*"1SD,- 


For  the  numerical  experiments  we  used  the  initial  conditions 


u(x,  0)  =  0 

4>(x,  0)  =  16  +  0.1  sin(27rx)  0  <  x  <  1 


and  took  =  16.  We  again  used  periodic  boundary  conditions  and  compared  Lax-WcndrolT  on  the 
unsplit  problem  with  k  —  h/20  to  the  split  scheme  with  k  =  h  on  the  slow  scale  and  the  method 
of  characteristics  for  Qf{k/2)  .  For  h  =  1/100  the  results  arc  shown  in  table  8.2.  Again  the  split 
scheme  outperforms  the  unsplit  scheme.  The  errors  were  reduced  by  a  factor  of  100  while  at  the 
same  time  the  work  was  reduced  by  roughly  a  factor  of  10. 
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0.5  C.S  1 


0.6  C.8 


0.6  0.8 


Figure  8.1.  True  and  computed  solutions  at  f  =  4.72  For  example  8.1.  The  first  component,  u\, 
is  on  the  left  and  the  second  component,  U2,  is  on  the  right.  The  schemes  used  are: 
top:  Unsplit  Lax-Wcndroff 
middle:  Time-split  method  (2.4a,b) 
bottom:  Time-split  method  with  change  of  variables 


Table  8.1.  Max-norm  errors  for  example  8.2  at  various  limes  t.  The  schemes  used  are: 
#1:  unsplit  Lax- Wend ro IT  with  k  =  h/2 
#2:  Leapfrog  Duhamcl  with  k  =  \h,  Q/{k)  =  exp(AraclI) 

#3:  Time-split  method  (2.1a, b)  with  k  =  ih 

#4:  Time-split  method  (2.4a,c)  with  k  =  ih,  m  =  8. 


6.6l9(-2) 
l.342(— 1) 
2.058(— 1) 
2.G85(—  l) 


l.336(— 2) 
l.949(— 3) 
1.414(— 2) 
3.434(— 3) 


2.147(— 3) 
4.598(— 3) 
7. 1 93( — 3) 
9.617(— 3) 


f>.470(-2) 
1 .3 15( —  l) 
2.0 1 6( —  1) 
2.623(— l) 


1/100 


1.677(— 2) 
3.389(— 2) 
5.314(— 2) 
6.071  ( —  1) 


3.356(— 3) 
4.130(— 4) 
3.365(— 3) 
2.028(— 4) 


5.58 1  ( — 4) 
1.166(— 3) 
l.845(— 3) 
2.437(— 3) 


1.635(— 2) 
3.320(— 2) 
5.197(— 2) 
6.818(— 2) 


Table  8.2.  Max-norm  errors  for  u  and  (/>  in  example  8.3  at  various  times  t.  The  schemes  used 
are: 

#1:  uiisplit  Lax-WendrolT  with  h  ■—  l/l 00,  k  =  h/20 
#2:  Time-split  method  (2.4a,b)  with  k  =  h  =  1/100. 


$ 


3.983(— 4)  2.952(— 6) 

3.354(— 5)  2.338(— 7) 


8.059(— 4)  5.882(— 6) 


0\ 
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